Pre- Turbulent Regimes in Graphene Flows 
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We provide numerical evidence that electronic pre-turbulent phenomena in graphene could be 
observed, under current experimental conditions, through detectable current fluctuations, echoing 
the detachment of vortices past localized micron-sized impurities. Vortex generation, due to micron- 
sized constriction, is also explored with special focus on the effects of relativist ic corrections to the 
normal Navier- Stokes equations. These corrections are found to cause a delay in the stability 
breakout of the fluid as well as a small shift in the vortex shedding frequency. Finally, a relation 
between the Strouhal number, a dimensionless measure of the vortex shedding frequency, and the 
Reynolds number is provided under conditions of interest for future experiments. 
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Since its recent discovery [HQ, graphene has continued 
to surprise scientists with an amazing series of spectacu- 
lar properties, such as ultra-high electrical conductivity, 
ultra- low viscosity to entropy ratio, combination of ex- 
ceptional structural strength and mechanical flexibility, 
and optical transparency. Many of these fascinating ef- 
fects are due to the fact that, consisting of literally a 
single carbon monolayer, graphene represents the first 
instance of a truly two-dimensional material (the "ulti- 
mate flat land" [3]). Moreover, due to the special sym- 
metries of the honeycomb lattice, electrons in graphene 
are shown to behave like an effective Dirac fluid of mass- 
less chiral quasi-particles propagating at a Fermi speed of 
about vf ~ c/300 ~ 10 6 m/s. This configures graphene 
as a very special, slow-relativistic electronic fluid, where 
many unexpected quantum-electrodynamic phenomena 
can take place, (j-[6|. In particular, the capability of 
reaching down viscosity to entropy ratios smaller than 
that of superfluid Helium at the lambda-point, has re- 
cently spawned the suggestion that electronic transport 
in graphene may support pre-turbulent phenomena, Ref. 

0. 

In this Letter, we pursue this suggestion in quantita- 
tive terms. More precisely, we simulate the relativistic 
graphene-fluid equations, proposed in Ref. [7], under 
conditions of present and prospective experimental re- 
adability. Our main result is that micro-scale impuri- 
ties, as small as a few microns, are capable of trigger- 
ing coherent patterns of vorticity in close qualitative and 
quantitative resemblance with classical two-dimensional 
turbulence (see e.g. Fig. []]). It is also shown that such 
vorticity patterns give rise to detectable current fluc- 
tuations across the sample, well in excess of flickering 
noise. As a result, based on our simulations, we con- 
clude that the hydrodynamic picture of graphene as a 
near-perfect, slow-relativistic fluid, as developed in Ref. 
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FIG. 1: Pre-turbulence at Reynolds number Re — 25 in 
graphene is shown: a t 37 940 tim e st eps, [(a)] and |(b)| and 
at 603400 time steps, [(c)] and [(d)] For [fb^] and [(d)] the term 
dp/dt was removed. The color represents the magnitude of 
the velocity. 



[7j, should be liable to experimental verification. The 
equations for the Dirac electron fluid in graphene read 
as follows [7]: dp c /dt + V • (p c u) = 0, for charge con- 
servation; de/dt + V • [(e + p)u\ = 0, for energy density 
conservation and 
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for momentum conservation. Here, c is the Fermi speed 
(~ 10 6 m/s), e the energy density, p the pressure, p c the 
charge density, and u the velocity. The shear viscosity 
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FIG. 2: Vortex shedding in graphene at Reynolds number 
Re = 100, using a grid of 1024 x 512 cells. The color scale 
represents the absolute velocity of the fluid. The picture was 
taken at 4 x 10 6 time steps. 



can be calculated by using [?J 
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where C v ~ 0(1) is a numerical coefficient, T is the tem- 
perature, a = e 2 /ehc is the effective fine structure con- 
stant, e being the electric charge of the electron, e the 
relative dielectric constant, and N the number of species 
of free massless Dirac particles. Additionally, the entropy 
density can be calculated according to the Gibbs-Duhem 
relation e + p — Ts. These equations have been derived 
under the assumption \u\ «c. 

The relativistic lattice Boltzmann (RLB), proposed by 
Mendoza et. al. [8|,|9j, is hereby adapted to reproduce, in 
the continuum limit, the equations for the Dirac electron 
fluid described above. The RLB model 18] was defined on 
a three-dimensional lattice with nineteen discrete veloc- 
ities. Since graphene is 2D, we have adapted the model 
to a two-dimensional cell with nine discrete velocities, 
linking each site to its four nearest-neighbors, four next- 
to-nearest neighbors (diagonal), plus a rest particle. Two 
distribution functions, fi and are used for the particle 
number and momentum-energy, respectively. These dis- 
tribution functions evolve according to the typical Boltz- 
mann equation in single-time relaxation approximation 
UM, fi(x + Sx,t + St) - fi(x,t) = -(/< - fD/r and 
gi(x + Sx,t + St) - gi(x,i) = - ^ q )/r, where r is 
the single relaxation time, and the equilibrium functions 
and g^ are defined in Ref. [8|, |9[ . The shear viscosity, 
according to this model is r] = (e + p) (r — \) c 2 St/3c 2 , 
where q = Sx/St is the ratio of the lattice spacing to 
time- step size. 

We choose the equation of state e = 3p, which de- 
pends on temperature in the relativistic regime, as e ^ 
(knT) 3 /(ch) 2 = T 3 (in normalized units c = h = fc# = 
e = 1) Thus, the shear viscosity r] would depend on 
the third power of the temperature, leading to a different 
relation than Eq. ((2]). However, in the Dirac fluid, the re- 
laxation time for the electrons depends on the inverse of 
the temperature, r re \ = (ha) 2 /ksT [12|, and, therefore, 
introducing this dependence into the relaxation time r of 
the numerical model, we obtain the correct function for 



the viscosity. In numerical units (Sx = St = c = 1), we 
set the relaxation time to r = tqTq/T + 1/2, where T is 
the initial temperature and To the initial relaxation time. 

The hydrodynamics equations are similar to the non- 
relativistic Navier-Stokes equations with the exception 
of the compressibility term ~ dp/dt. This term is most 
likely negligible at low frequencies, but it may become 
relevant at higher ones. The Reynolds number Re, mea- 
suring the strength of inertial versus dissipative terms [?J , 
is given by Re = (sT/c 2 )(LUt yp /r]), where L and Ut yp are 
the characteristic length and flow velocity of the system, 
respectively. In lattice units, it reads as 
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According to classical turbulence theory, vortex shedding 
in graphene is expected for Reynolds numbers well above 
one, typically Re ~ 10 100. To detect signatures of 
pre- turbulent behavior in graphene experiments, one can 
measure the fluctuations of the electric current through 
the graphene sample. The current density is defined by 
j = p c u, and the total electric current is calculated inte- 
grating the current density along the transverse (y) co- 
ordinate. The characteristic fluctuation frequency can 
then be related to the vortex shedding frequency. Macro- 
scopic speeds u ~ 10 5 m/s could be achieved by the elec- 
trons in graphene [l3| . The Reynolds number rewrites as: 
Re = UtypLT j c 2 (rj / s) . According to Ref. [7], r]/s takes 
values around 0.2/i/fcs, at temperature of 300K, so that 
we can write Re = Ut yp L/v e ff, v e ff = c 2 r]/Ts ~ 0.005 
m 2 /s being the effective kinematic viscosity. Therefore, 
a sample of size L = 5/im, within reach of current tech- 
nology, would yield Re ~ 100, sufficiently high to trigger 
pre-turbulent phenomena, such as vortex shedding. To 
test the idea on quantitative grounds, we implement a 
simulation on a grid with 1024 x 512 cells. The follow- 
ing initial values (numerical units) were used: e = 0.75, 
p c = 1.0, u = (u x ,0) = (0.002,0), and the Fermi speed 
c = q = 1.0. The initial value of the relaxation time was 
chosen To = 0.003 such that the initial shear viscosity 

, = H--^) = io- 3 . 

A circular obstacle, with diameter D = 50, is intro- 
duced at (256, 256), modeling a 5 micron diameter impu- 
rity in the graphene sample (Fig. [2}. With this configu- 
ration, and setting L = D in Eq. [3j the Reynolds number 
for this system is Re ~ 100. We choose periodic bound- 
ary conditions at top and bottom, and demand that the 
distribution functions of the boundary cells are always 
equal to the equilibrium distribution functions evaluated 
with the initial conditions. Free boundary conditions are 
imposed at the outlet. At the left border, we set inlet 
conditions, where the missing information of the distri- 
bution functions is filled by the equilibrium distribution 
function corresponding to the initial conditions [14]. We 
define £7; = 0.05ps. 
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FIG. 3: FFT of the electric current fluctuations AI g , in the 
graphene sample (top right), due to the vortex shedding as 
a function of time. Also, it is shown the FFT of the drag 
force acting on the obstacle (top left). This result refers to 
Re = 100. At the bottom, the fluctuations in the electric 
current I g are shown as a function of time. 



The drag Fjj x and lift Fjj y forces acting on the obsta- 
cle are measured, the vortex shedding frequency being 
computed in terms of fluctuations of the lift forces. We 
compare the frequency of the electric current fluctuations 
with the frequency of the drag force, which, in general, 
is twice the vortex shedding frequency (see Fig. [3j). To 
relate these to the vortex shedding, we use a fast-Fourier 
transform (FFT). As is well visible from Fig. [3j the cur- 
rent fluctuations contribute about one part per thousand 
of the base signal, and, consequently, they should be li- 
able to experimental detection. In future applications, 
involving larger graphene samples, higher Reynolds num- 
bers will be attained. Consequently, it becomes of inter- 
est to assess the role of the relativistic corrections to the 
classical Navier- Stokes equations. 

Comparing the dynamics of the relativistic and non- 
relativist ic fluids, two basic differences emerge: the rela- 
tivistic correction term ~ dp/dt; and the viscosity depen- 
dence with the temperature, Eq. (|2j). In order to assess 
whether these terms play an important role, we imple- 
ment three simulations on a grid of size 2048 x 1024 cells. 
In the first simulation, we model the full relativistic equa- 
tions; in the second one, the relativistic effect ~ dp/dt is 
removed; and in the third one, the viscosity is forced to 
be a constant. The same initial configuration, as before, 
is used with the exception of: u = (u x , 0) = (0.03, 0) and 
D = 100 (in this case, modeling an impurity of diame- 
ter 150/im). The impurity is now centered at (512,512). 
With this configuration, Eq. [3] gives Re ~ 3000. The 
simulations run up to 10 6 time steps (with 5t = 1.5ps). 

From Fig. [H we find that, in the case of constant vis- 
cosity, the frequency is a bit higher than the one corre- 
sponding to the full relativistic case. On other hand, if 
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FIG. 4: Frequencies of the vortex shedding at Reynolds num- 
ber Re = 3000, using a grid of 2048 x 1024 cells, are shown. 
These are calculated for three different cases: the full rela- 
tivistic, relativistic without the term dp/dt, and relativistic 
with constant viscosity. 
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FIG. 5: Strouhal number St as a function of the Reynolds 
number Re for both non-relativistic and relativistic fluids. In 
the inset, the mean value of the x-component of the drag force 
as a function of the Reynolds number is shown. The error bar 
is of the size of the symbol. 



the term ~ dp/dt is removed from the equations, the fre- 
quency decreases. We conclude that, in order to compare 
to high precision measurements of the vortex shedding 
frequencies, these terms cannot be ignored. To study 
the frequency of the vortex shedding, we vary the initial 
velocity in order to obtain different Reynolds numbers. 
The Strouhal number St is defined as the dimensionless 
frequency of the vortex shedding and can be calculated 
as St = f s L/Utyp, where f s is the frequency of the vor- 
tex shedding. Fig. [5] shows that the relation between 
St and Re is very similar for the relativistic and non- 
relativistic fluids, with a fast growth of St in the range 
200 < Re < 1000, followed by a flat-top at St ~ 0.25 
151-17] for Re > 1000. From the Strouhal number, we 
can obtain the frequencies of the vortex shedding, as 
f s = 0.2Ut yp / L. The frequency of the drag force is twice 
that of vortex shedding, namely f^s — 0AUt yp / L. As 
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FIG. 6: The same case of Fig. [4) but in the case of the con- 
striction at Reynolds number Re = 25, using a grid with 
1024 x 1024 cells. 




~ 1 2 4 6 8 10 12 

Time (ns) 

FIG. 7: The same case as Fig. [3] for the constricted flow at 
Re 25. 
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FIG. 8: The same as Fig. [5] for the case of the constriction. 



a result, once the Reynolds number is known, one can 
compute the frequency of the drag force, the Strouhal 
number, and then compare with the FFT of the electric 
current measurement in the sample. The mean value of 
the drag force Fd x , reported in the inset of Fig. [5] as 
a function of the Reynolds number, shows a monotonic 
dependence in the range of Re explored here. 

Another kind of set-up to detect pre-turbulence in 
graphene experiments, with the possibility of being im- 
plemented nowadays, consists of building a constriction, 
where the Dirac fluid can develop vorticity patterns as it 
crosses through. Fig. [T] shows the vorticity at Re = 25, 
where the characteristic length L = 50 cells has been 
chosen as the distance between the tips. In this case, 
the initial velocity is taken u = (u x ,0) = (0.0005,0), 
in lattice units, and the simulation is performed using 
a grid of 1024 x 1024 cells. We simulate two systems, 
one with the full relativistic equations and the other one 
by just removing the relativistic term dp/dt. From the 
simulations (see Fig. [I]), we conclude that the relativistic 
contribution affects the time to the onset of instability, 
and, from Fig. [6j we can appreciate that, as for the cir- 
cular impurity, the frequency of the vortices presents a 
shift due to the relativistic corrections. However, both 
constant viscosity and removal of the relativistic correc- 
tion, contribute to an increase of the frequency of the 
fluctuations. Fig. [7] shows how such fluctuations can be 
measured, and the characteristic frequencies (see red cir- 
cles in Fig. [7j) related with the drag force acting on the 
constriction. Note that, in order to achieve Re = 25, at a 
speed of 0.1c, the distance between tips is about 1.25/im. 

As for the case of the circular impurity, we can find 
the characteristic relation between the Strouhal number 
and the Reynolds number for this geometrical set-up (see 
Fig. [8]). From the inset of Fig. [5J we observe that the 
drag force decreases slightly, as the Reynolds number is 
increased, and exhibits a noticeable difference between 
the non-relativistic and relativistic cases. 

Summarizing, we have shown that, in the range of 
Re ~ 10 2 , vorticity patterns can be indirectly observed 
by measuring the electric current fluctuations in the 
graphene sample. However, using a different geome- 
try, like a constriction, signatures of pre-turbulence can 
be detected already at Reynolds numbers as small as 
Re ~ 25. We have also compared the effects of relativis- 
tic corrections, such as dynamic compressibility and the 
dependency of the viscosity on the temperature, on the 
dynamics of the system. In these cases, the temperature 
dependency of the viscosity and the term dp/dt produce 
a shift in the frequencies of the vortex shedding and, 
therefore, in the electric current fluctuations. Addition- 
ally, the relativistic correction term, ~ dp/dt, is found 
to delay the instability process in the case of the con- 
stricted flow. For future applications, most likely access- 
ing higher Reynolds numbers, the frequency of the vortex 
shedding can be calculated using the Strouhal number, 
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thereby permitting to distinguish current fluctuations in- 
duced by pre-turbulent phenomena from those resulting 
from other physical effects. 

We thank K. Ensslin for enlightening discussions. 
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